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ABSTRACT 


The flow field and the temperature distribution around a heated/cooled 
circular cylinder placed in an insulated vertical channel is determined using 
a novel finite volume algorithm which combines some important features of 
the finite element methods with the finite difference based techniques. Ele- 
mentwise interpolation (iso-parametric) and transformation of non-orthogonal 
element geometry into a square computational element have been made use 
of while solving the integral conservation equations. Unsteady mixed convec- 
tion situation has been considered for the study and the effect of buoyancy on 
the vortex shedding in the wake of a cylinder has been observed for a fixed 
Reynolds number of 100. Vortex shedding at Re=100 is found to stop com- 
pletely at a critical Richardson number of 0.15. Below this critical value of 
Richardson number, shedding of vortices into the stream is quite prominent. 
Heat transfer characteristics of the heated/cooled cylinder are studied as a 
function of Richardson number. The results are of use in the design of shell 
and tube heat exchangers. 
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Chapter 1 

INTRODUCTION 


1.1 Flow past bluff bodies 

The flow past bluff bodies has been a fascinating problem both for aca- 
demic and industrial applications. The behaviour of the flow in the wake of 
a bluff body is responsible for changes in heat transfer and flow induced vi- 
brations on the body. Flow separation and vortex shedding tend to increase 
the heat transfer but at the same time increase the drag on the body. The 
shedding may prove to be detrimental for the life of the body. The objective 
of getting maximum heat transfer without being affected much by the flow 
induced vibrations provides with enough challenge in this field of study. 

Understanding of the flow past a bluff body is useful for the design of heat ex- 
changers, air conditioning coils, nuclear reactor fuel rods, chimney stacks, T.V 
towers, suspension bridges etc. The Reynolds number associated with these 
flows are sufficiently high to produce turbulent flows in the wake region even if 
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the boundary layer remains laminar. Slaouti and Gerrard(1992) have pointed 
out that there is no general turbulence model to incorporate into the numerical 
simulation of unsteady separated flows. Direct numerical simulation of Navier- 
Stokes equations is possible but it requires large computing power. It has been 
observed in experiments that some features of large scale wake behaviour as- 
sociated with vortex shedding at high Reynolds number are similar to that 
at low Reynolds numbers. At low Reynolds numbers, predicting flows and 
heat transfer around a cylinder may be reasonably good by two-dimensional 
laminar computations. 

1.2 Mixed Convection 

The diiTerencc in temperature of the bluff body and the fluid in many 
practical situations demand the study of mixed convection. In a vertical flow, 
the velocity induced by free convection will be either aligned or opposed to that 
of the forced convection. It hcis been observed that in the presence of aligned 
free and forced convection, the wake region is diminished and the separation 
points are shifted downstream. When the free convection acts in the opposite 
direction the point of separation is advanced upstream. 
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1.3 Literature Review 

Fromm and Harlow [1], Thoman and Szewczyk [2] and Jain and Goel [3] have 
studied numerically flow past a circular cylinder. Forced convection over a row 
of inline cylinders placed between parallel plates has been studied by Kundu et 
al [A, 5]. The effect of buoyancy on the flow structure has not been taken into 
consideration in the above investigations. When flow velocity is not very high 
and the temperature difference between the bluff body and the fluid is signifi- 
cantly high, the heat transfer characteristics are strongly influenced by thermal 
buoyancy effects. In many practical applications, the cylinder is placed inside 
a confined space so that the side walls have considerable influence on the heat 
transfer rate. Quite a few efforts dealing with mixed convection regime have 
used steady state equations and arc limited to small Reynolds numbers up to 
40 (Farouk and Guceri [6] and Badr [7]). 

Biswas et al [8] have studied unsteady mixed convection heat transfer in a 
horizontal channel with a built-in square cylinder. They have found that for , 
a given Reynolds number, the heating of the fluid in the channel is improved 
by mixed convection up to a certain Grashoff number and deteriorates if the 
Grashoff number is further increased. Oosthuizen and Madan [9] studied ex- 
perimentally mixed convection for Reynolds number ranging between 100 and 
300. Noto and Matsumoto [10, 12] reported degeneration of unsteady vortices 
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to the stationary ones and termed it eis the breakdown of von-Karman street. 
Chang and Sa [13] have studied the behaviour of the near wake vortices. They 
have predicted the degeneration of purely periodic flows into a steady vortex 
pattern at a critical GrashofF number of 1500. 

The present study considers an unsteady mixed convection flow over a cir- 
cular cylinder placed inside an insulated vertical channel. In order to compare 
the results of present investigation with that of Badr [7] and Dennis et al [19], 
some computations have been performed for the situation concerning flow over 
a heated cylinder placed in an infinite medium. Extra-Flag algorithm proposed 
by Mukhopadhayay et al [14] has been used for the analysis. 




Chapter 2 


STATEMENT OF THE 
PROBLEM 

2.1 Domain Definition 

The system of interest consists of a vertical channel with insulated side walls. 
An obstacle in the form of a circular cylinder is placed inside it (Fig.l). The 
blockage ratio (D/I I) is equal to 0.25 and the length of the channel is twenty 
four times the diameter (D) of the cylinder. The centre of the cylinder is placed 
on the vertical axis of the channel at a downstream distance of 8D from the 
entry plane . 

2.2 Governing Equations 


The domain is descretized into small curvilinear quadrilateral cells (Figs. 2, 3). 
The velocity and temperature nodes are located at the vertices, and the pres- 
sure nodes are located at the centroids of these cells. The continuity control 





(a) Complete channel 


(b) Near the cylinder 


Fig.2 Grid 




CVl 

mass continuity 
control volume 

CV2 

momentum 
control volume 

• - Velocity and temperature 
X -Pressure 

(a) Schematic representation of tlie discretization 


Nk Uk 


NkVk 

Mk Pk 
k=l 



Interpolation domain 
for velocity 


Pressure node 


Central velocity 
node 


* 

Integration domain (CVS) 



(b) Interpolation domain for momentum balance 


Fig. 3 Domain discretization 
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volumes (CVl) are formed by the grid lines connecting velocity nodes and the 
momentum control volumes (CV2) are formed by curvilinear quadrilaterals 
whose vertices are the pressure nodes . 

Considering a typical control volume for mass continuity, CVl, the integral 
form of the mass balance can be written as 

I- / / pdA+ f pV-ndl = 0 (2.1) 

Ot J Jcvl Jol 

where, dl is an elemental length on the boundary CSl of the control volume 
CVl, and n is the local outward normal vector. For incompressible flow, the 
first term of the equation(2.1) is identically zero, which simplifies the equation 
to 

f V.' ndl = f [urix + vny)dl = 0 (2.2) 

Jcsl Jcsl 

where, rix and Uy are the direction cosines of the outward normal n on the 
boundary CSl. 

Now considering any typical momentum control volume CV2, the integral 
form of X and y momentum balance equations can be written as 

■^ [ f u{pdA) + f u{pV-hdl) = Fx (2.3) 

Ot J Jcv2 J ca2 

^ff v{pdA)+ f v(pY_-hdl) = F, (2.4) 

Ot J Jcv2 Jca2 

where Fx and Fy are the components of the resultant forces acting on the 
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control volume in x and y directions respectively. These can be evaluated as 


F. 


= i • I a - ndl 
J cs2 

= t • j^^{-pi+Ai(vK + yKl)- 


(2.5) 


where a is the stress. For incompressible flows 


F. = J (- 

Jca2 I 


^ du f du dv 




dl 


( 2 . 6 ) 


For the y momentum equation, 


Fy = j • f g:-hdl +j ■ [ f pgfi{T - Too)dA (2.7) 

J C32 j j cv2 

where T is the absolute temperature. The second term in the above equa- 
tion is due to buoyancy force. For incompressible flows 


F = 

X 7/ — 


f ( n 1 Jt 

L[-Pn, + + p 

// pgP(T-T^)dA 

J J cv2 


+ 


( 2 . 8 ) 


Finally the x and y momentum equations can be written as 


^ 2 ^ I 2 vny)(ii| = 

L r'”*' + + " (a? + & j J 


(2.9) 
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" {I JL + L + ''"-'>‘"1 = 

f ( n 1 11 

Jcs2 i ^ ^ j 

+ // pg(3{T-T^)dA (2.10) 

^ ^/ct;2 

Now considering a typical control volume, CV2, the integral form of the 
energy equation can be written as 

o 

<pdA) + ^(Ze) • ndZ = ^ kVT • Ml + J ^ QdS (2.11) 

where e is the internal energy, k is the thermal conductivity , T is the 
temperature of the fluid, and Q is the volumetric rate of heat generation. For 
the case with no heat source, Q is identically equal to zero. For incompressible 
flows. 


cT = e 


where c is the specific heat and it is a constant (which is a reasonable assump- 
tion for incompressible flows). With- this substitution, the energy equation' 
becomes 


[I L + L = L ^ f "“} 
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The non-dimensionalized form of the equations is as follows: 


/ {Un^ + Vny)dl = 0 

Jcsl 


(2.13) 


dr 


J I + I (Un. + Vny) dlj = 


Ic32 { 


„ 2 dU 1 (du dV\ , 

Ptlx “hr. o Tr “h r» lo-ir"h — 1 ri„ / dl 


RedX " Re \dY dX 




(2.14) 


dr 


If VdA+ I V{Un, + Vny)dl\ = 

J J cv 2 J cs 2 I 


<fc32 ^ 


p ^ 2 av 1 (du ^dv\ \.. 

Re ^ Re dX ) J 

Gr 


n 2 f f 

Re^ J Jcv2 


(2.15) 


dr 


J I + I + Vny) d/| = 

I f 1 

Jc32\PrRe \dX^" 


39 \ 

+ ^ny dl 


(2.16) 
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2.3 Boundary Conditions 

Parabolic velocity profile at ambient temperature is assumed at the inlet. The 
channel walls are insulated. The cylinder temperature is governed by the 
magnitude of the Richardson number. No-slip boundary condition is imposed 
both on the cylinder and the channel walls. The exit boundary is located at far 
downstream so that it does not influence the vortices in the near wake region. 
However, at the outlet, the derivative of the variables is made equal to zero. 
Mathematically one can write 

at the inlet 


U = 0 = 0-, V= V{X) 


at the outlet 


dY~^' ar ~ ° ’ dY~^ 


at the channel walls 


U = V = 



at the cylinder surface 
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[/ = 1/ = 0 ; 0=1 

Some cases have been computed to simulate the flow past a circular cylinder 
placed in an infinite medium. For such a simulation, the blockage ratio (D/H) 
has been taken as 0.1. The boundary conditions at the side walls have been 
modified to simulate the far field boundary. The boundary conditions for the 
infinite medium situation are as follows ; 
at the inlet 

U = 9 = 0 ; V = 1.0 

at the outlet 

^-n —-n 
dY~^' dY~^' dY~^ 

on the sides 

C/ = <? = 0 ; y = 1.0 

at the cylinder surface 


U = V = 0] 0 = 1 




Chapter 3 

METHOD OF SOLUTION 

3.1 Description of the Extra-Flag algorithm 

The EXTRA-FLAG algorithm proposed by Mukhopadhyay et al [H] has 
been extended for the present work. The energy equation is solved together 
with the Navier-Stokes equation. EXTRA-FLAG algorithm is based on in- 
tegral maiss, momentum and energy balance applied over the non-orthgonal 
control volume. The integral flow equations are converted into algebraic form 
through numerical quadrature coupled with spatial interpolation (isoparamet- 
ric) of their nodal values. In order to facilitate the numerical quadrature, each 
non-orthogonal quadrilateral control volume is mapped into a standard square 
cell with the help of the local curvilinear coordinates. Explicit time marching 
is employed for the momentum equations using the pressure field correspond- 
ing to the previous time level. Principle of elemental mass conservation, on 
the other hand, is enforced in an implicit manner by iteratively correcting the 
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velocity and pressure fields. The method of solution thus utilises the geomet- 
rically additive features of FEM ( isoparametric interpolation and numerical 
quadrature) with the control volume approach applied in MAC [15] and SIM- 
PLE [16] family of algorithms. 

Since the order of the momentum equations is one order higher than the 
continuity equation, the velocities have been interpolated using bi-quadratic 
functions, while applying the momentum principle. For the satisfaction of con- 
tinuity, bilinear interpolation functions have been used to describe the velocity 
variables over the respective control volumes. The temperature interpolation 
is similar to that of the velocities for the momentum control volume. It can 
be represented mathematically as follows : 

For continuity equation 


U = j2MkUk-, V = (3.1) 

it=i fc=i 

For momentum and energy equation 


U = j2NkUk; V = f2NkVk; e = j^Nkek (3-2) 

fc=i *=i fc=i 

The summation in equation (3.1) is performed by mahing use of velocities 
at the four vertices of the continuity control volume. Equation (3.2) uses the 
velocities of the neighbouring nodes of the corresponding momentum control 
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volume. These are shown in Fig.3. The variation of the pressure over the 
momentum control volume is expressed in the form 


P = f^MkPk (3.3) 

k=i 

In the equations , Mk and Nk are the interpolation functions, and are given 


by 


w. = - m - 2) 

JV, = - 2)(( + 2) 

= ^<’)(’I + 2)K+2) 

ffr = + 2)(( - 2) 

JV, = ^(4 - r,^)(4 - e) 
M: = ^(1 - 0(1 - 1) 

M , = ^(1 + 0(1 + r ,) 


N2 = - e)(v - 2) 

Nt=y(4-ri0(( + 2) 
N, = T,(4 - e)(r, + 2) 
Wb = ^{(4 - >)")({ - 2) 

and 

M, = i(l + {)(! - v) 

M, = i(l - fl(l + r,) 


( 3 . 4 ) 


The integrals appearing in equations (2.13)-(2.16) have been evaluated us- ' 
ing Gauss- Legendre quadrature. This is done by mapping the non-orthogonal 
control volume into a square computational cell whose co-ordinates vary from 
-1 to 1 in both the directions (Fig.4). The mapping is done by isoparametric 


transformation. 




Fig.4 Transformation of a curvilinear control volume 
into a square computational cell. 
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The co-ordinates x and y for any point within the non-orthogonal control 
volume are interpolated in terms of their nodal values in much the same way 
cLS the flow variables. The interpolation scheme for the control volume CV 2 
arc given by 

9 9 

X = '£ NkXk ; F = X] NkYk (3.5) 

*=1 jk=i 

where {Xk,Yk) are the co-ordinates of the nodes surrounding the control 
volume as shown in Fig 3(b). For continuity cells (CVl), the linear transfor- 
mation is used in the form 


X = J2 ^kXk ; F = E MkYk (3.6) 

k=:l k=l 

With the above transformation the integrals appearing in the equations 
(2.13)-(2.16) can be rewritten in the form of transformed co-ordinates ^ and 77 . 
For all the mapped cells, ( and 77 vary from -1 to 1 , and the Gauss- Legendre 
quadrature can be easily adopted for evaluating the integrals. The area inte- 
gral is given by 


/ f Fa,ri)d(dr, = (3-7) 

7 - 1 7-1 t=, 

where {^kiVi) ^^e the Gauss points and Wk, and Wi are the corresponding 
weights. The co-ordinates of the Gauss sampling points are the zeroes of the 
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degree Legendre polynomials, where m is the total number of Gauss points 
used in a particular direction. In the present case it ie chosen to be equal to 
2. The line integral can be evaluated as 


rl 

/ nm = E 

k=l 

rl ^ 

/ F{r})dT] = X^F(77fcVfc 


For the completion of the transformation, it is neccessary to convert the arc 
lengths, areas, surface integrals, derivatives and integrals in the physical co- 
ordinate system to the computational co-ordinates ^ and rj. For the typical 
control volume CV2 (Fig. 4) chain rule may be applied in the following 


X = Xi^,rj) ; Y = Y{(,rj) 


dX dX j 

dx = ^di^~d-n ■■ dY 




or in the matrix form 


■ dX ■ 


1 






. dY _ 


1 

1 


drj 



l^k^l J 


. 


in a similar fashion, the derivatives can be obtained as 
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r dx 

dY 1 




r d 1 



Bi 


r ^ 1 


r ^ 1 

dc 





dX 

= [J] 

dX 

d 





d 


d 

L dv - 


dX 



L dv - 


1 dY J 


- dr] 

dr) _ 





where [J] is the Jacobian matrix for the transformation. The matrix in 
the above equation can be evaluated using equations (3.4) and (3.5) in the 
definitions given above. From the equation (3.8), the velocity derivatives can 
be written as 


dU 1 


r dU -1 


r Ml 

dNi 

dN^ 1 


\Ur] 

dX 

= [J]-' 

B( 

= [jr 

Bi 

Bi 



U 2 

dU 


§!l 


dNi 


m. 


1 

dY J 


- dr) - 


- dr) 

dr) 

dri - 


1 


Similarly we can evaluate the derivatives of V,(?. 


On the sth side (s=l,2,3,4) the elemental length vector dl, is 


du 


dX,i -b dY ,j 



(3.10) 
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The outward unit normal vectors on the corresponding cell boundaries are 
given by 


n, — 






fa] 

\9CJi 


(3.11) 


where m = 0 for s=l,2 and m=l for s=3,4 and 


■ \d< 




(3.12) 


C = ^ for s=l,3 and ( = rj for 8=2,4. The elemental length, dl is repre- 
sented zis 


dl = di (3.13) 

for the sides 1 and 3, since rj is constant for these sides. Similarly 

= (^) 

for sides 2 and 4, since ^ is constant for these sides. 

Lastly , the elemental area dA can be evaluated as 

dA = \J\d^dT] (3.15) 

For continuity control volume, it is not necessary to employ the numerical 
quadrature, since the shape of CVl is simpler and the velocity interpolation 
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is linear. 

Using the interpolation as described above the final form of momentum 
and energy equations can be written in the matrix form as 


where 


f ^ 1 

1 

f U ) 


U[[C] + [D]] V 

1 e J 

1 

1 ^ J 



o 

o 


[C] = 

0 C2j 



0 0 


-t- [S] [P] -b [Sr] 


0 


"3j J 


and the elements Cij, C 2 j and Caj are given by 


(3.16) 


(3.17) 


[^ij] = [^ 2 j] — [Caj] — -z f Nj{Unx + Vny)dl — 


-1 


4 2 

EE 


cv2 

9 


Nj X^(t7,n^ -f- Vgny)Ng 
9=1 


The matrix [D] in equation (3.16) is given by 


Wk^ls 


[D] 


DUjj DVij 0 
0 

0 0 DUaj 


DUaj DVaj 


(3.18) 


(3.19) 


where 
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DUu = 


I r ( 2 dNi . .... . 

lcv2 L2 \Re dX Re dY 


DU2j = 


DVu 


DO: 


3i 


1 dNj 


■EE 


2 dNj 1 a/Vj 


Ac „2 t'l aA: ^ Re dY 


n. 


WkAlg 


(- 

^C‘u2 cs2 y. 


1 dN 


frij; dl 


-n. 


= 


A(;^2 V/?e ay 
1 ^ ^ r 1 aiVj 
Acv 2 hi hi ay 
J_ f fj^9Nj 

Acv2 Jcs2\ 

Act/2 

A 


u^fcAia 




i^eax”*' 


EE 

s=i jk=i L 

^ r fj^dNj 

.cv2 Jcs2 \ 


WkAlg 

2 dNj 


n„ dl 


1 


A 


EE 

c^'2 s=i k=l 


Re ax ^ ReW''\ 
2 r 1 aXj-^ 2 dNj 

Re dx ay 


lojfeA^, 


-f —( 

lct;2 -/ciZ Pr/le \ 


ax,’ ^ax,- \,, 

M"-+sF"»r' 


EE 


Ac ,;2 PrRe 


dNj dNj 

+ -TrA-Ur 


dX 


dY 


WkAl, 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


(3.24) 


The pressure matrix in turn is 


where 



(3.25) 


Sij = —j— / = —T~ E 

^cv2 ^ ^cv2 5=1 jt— 1 

S2j = —h~ / = -7^ E 

Act;2 -^cv2 ,=1 fc=l 


(3.27) 
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I^astly, the buoancy term [Sr] can be expressed as 


where 


[Sr] = 


0 

Sra 

0 


(3.28) 


II ( 3 - 29 ) 

ylcv2 -I Jcv2 lic^ 

where O^v is the non-dimensional temperature at the center of the control 
volume CV2. 

The continuity equation can be discretised as 


-h = 0 (3.30) 

5=1 

3.2 Upwinding of Convection Terms 

In the numerical solution of Navier-Stokes equation, the convection terms need 
specific attention for conserving the transportive property. An effective upwind 
discretisation for the convection terms is implemented by setting 


C.jUj = C 2 jUj = -^ I {Un, + Vny)GjUjdl (3.31) 

^cv2 

where Gj is the upwind basis function, and is given by 


central jgRAB' 


P'^ n.lVAT ION O F PR£-SS_U RE CORRECTION EQUATION M 

E GiVi = E(1 - a)NjUj + aU, (3.32) 

The subscript b corresponds to the node situated upstream of the concerned 
face and a. is the upwinding factor. For purely symmetric interpolation a = 

0 and for the upwind biased interpolation, a = 1. In all my computations, I 
have used o;=0 so that there is no artificial viscosity. 

3.3 Derivation of Pressure Correction Equa- 
tion 

When the time gradient of velocities and temperature arc expressed in the 
discretized form, equation (20) becomes 

I }+A^1C” + D”)| | + Ar[S){P”+'} + |Srr (3.33) 

where the subscripts n and n + 1 denote time levels. However, at the start of 
calculation for each time step, is not known and hence equation (3.33) 
predicts provisional velocities (represented by a superscript asterisk *) that 
satisfy momentum balance corresponding to the pressures. When both 
pressure and velocities undergo iterative corrections, F" takes up provisional 
value P*. Thus 


{ ((; } = { ((I } + At 1C” + D“] { (^! } + At (S] {P'} + [S]" (3.34) 
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As the continuity equation is satisfied in each cell, all provisional velocities 
reach the final values for the time step (n+1). Subtracting equation (3.34) from 
equation (3.33) the expressions for the velocity corrections can be obtained as 

where 5U, SV and 5P are the velocity and pressure corrections, respectively. 
Applying the principle of continuity at the (n + l)th level, 

E M, = 0 (3.36) 

S=1 

The unconverged velocities, however, do not satisfy the continuity equation, 
leading to a non-zero residue of the form 

• E ^ 

5 = 1 

Subtracting equation (3.36) from equation (3.37), we obtain 

E 

5 = 1 

where SUg and 6 V 3 are the average corrections on the sth side, so as to sat- 
isfy mass continuity in the control volume. The quantities SUa and SV, are 
expressed in terms of the nodal velocity corrections, which in turn are substi- 
tuted by proper pressure corrections using equation (3.35). These average 
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velocity corrections, when evaluated as the simple average of the nodal 
velocities, lead to decoupling of pressure in alternate continuity cells. This 
decoupling of pressure variable was also observed by Majumdar et al [17]. 
Introducing face-center velocity corrections [//,,, V},, as shown in Fig.5, the 
average corrections on the east side of the CVl can be expressed as 


5U, = 

hc/„ + iw,„ + Uu,., 

6 6 6 

(3.39) 

= 

-SVll -f jSVllJ ■+ -iJV/.e 

6 6 6 

(3.40) 


where II and III are the velocity nodes lying on the east face. Similarly, 
the average corrections on the other sides can also be estimated. Note that 
equations (3.39) and (3.40) amounts to the application of Simpson’s rule along 
the boundaries of CVl. In order to find the expression for a 

modified form of the equation (3.35) can be written for the velocity correction 
vector as 


5V = -5tA{5P) 

eventually we can write 


(3.41)' 


SU 


I.'- 


(i . SV)j,^ = -lri- {8 ’J^} 




sp e. - IPe\ 
Me i 


(3.42) 



3.3. DERIVATION OF PRESSURE CORRECTION EQUATION 


30 


Similarly, 





(3.43) 


Here e is the unit vector connecting the central pressure node P to the eastern 
neighbouring pressure node E. Similarly, the face-center velocity corrections 
can be obtained as 


6U 


f,(n,w,s) 


At 


^(N,W,S) — \ — SP{N^w,S) \ 


^1{N,W,S) 


Al. 


(N,W,S) 


7 


I Ai(Ar,w,5) J ^ Af(iv,VK,5) J 


(3.44) 


where the lower-case subscripts (n,io,s) refer to the corresponding boundary 
faces, while the upper-case subscripts {N, W, S) refer to the corresponding 
continuity cells. The mass residue contribution to the east face is given by 


Ale{SUen^,c + ^Kny,e) = -Re 


(3.45) 


Substituting for 5Ue and 5Ve from equations (3.39) and (3.40), we have 
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A4 + ^iOui + n,,,+Al, + giV),} n,,. = -R, 

(3.46) 

From equation (3.35), the nodal velocity corrections are given by 

I I = At[S]„ {ips SpsE Sps SPF f (3.47) 

I I = Ar[S]ifj SpfiE SpB ^Pp (3.48) 

It is to be noted that in equations (3.47) and (3.48) the pressure coefficient ma- 
trix S is evaluated at the nodes. After substituting equations (3.42), (3. 43), (3.47) 
and (3.48) into equation (3.46), the residue contribution of the ea»st face, Re, 
can be calculated in terms of pressure corrections. Similar exercise for the 
other sides of the continuity cells would lead to a pressure correction equation 
of the form 


[CP] {SP} = -R (3.49) ’ 

« 

where CP is the coefficient matrix corresponding to the pressure correction 
array and R is the mass residue arising from velocities updated by momentum 


equations. 



Chapter 4 

RESULTS AND DISCUSSION 


Air (Pr=0.7) has been used as the working fluid. Boussinesq approximation has 
been applied while writing the momentum equations The results are obtained 
for a Reynolds number of 100 and for' different values of Richardson number. 
Streamlines for various Richardson numbers have been plotted in Fig. 6. These 
plots correspond to a dynamic steady state behaviour. The physical signifi- 
cance of negative Richardson number is culmination of a buoyancy force in a 
direction opposite to the flow direction. In the viscous region, at the vicin- 
ity of the aft of the cylinder, the inertia force is opposed by the force due to 
adverse pressure gradient, the negative buoyancy force and the viscous force. 
Depending on the magnitude of Richardson number, the point of separation 
moves towards the leading edge causing an early separation. This also broad- 
ens the wake. Vortex shedding becomes quite prominent in the entire regime 
of negative Richardson number. The shedding phenomenon continues 
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to c.haractcrizo the; flow even at lli=0 and Ili=0.1. Fur a Ridiardson num- 
ber of 0.15(not shown here) and above, the flow structure becomes different, 
'riu; aiding buoyancy changes the dynamics altogether. Near the vicinity of 
the cylinder, even at the aft, the buoyancy force is added up with the inertia 
force in the face of pressure force and viscous force. I’his results in a sepa- 
ration delay and the vortex shedding is also stopped. This observation is in 
confirmation with the earlier studies by Noto and Matsumoto [10, 11, 12] and 
(Jhang and Sa [13]. h'igure 7 shows isotherms for various llichardson num- 
bers. The buoyancy opposed and buoyancy aiding cases are discerned in this 
figure. It can readily be appreciated that the waviness of the isotherms in the 
downstream of the cylinder is more for the buoyancy opposed convection as 
compared to buoyancy aiding convection. Beyond a Richardson number of 0.2 
the plume spread becomes narrow and waviness disappears. Figure 8 shows 
the phase diagram for different Richardson numbers. For all the Richardson 
numbers, the transverse velocity component at a point 3D downstream of the 
cylinder has been plotted against the streamwise velocity component at the 
same point. These plots clearly establish the periodic nature of the flow for all 
the cases below a Richardson number of 0.15. The orbits of the phase-diagram 
also indicate superimposition of oscillations with two different frequencies. 
The signal traces of the transverse velocity components and their FFl for the 
Richardson numbers of -1.0 and 0 have been shown in Fig 9. For the case 
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of Ri — -1.0, the peaks at St — 0.2861 and St=0.8583 confirms the presence 
of multiple frequencies (super harmonics). In the case of Ri=0,'we observe a 
similar trend which can he depicted as two component tiioe series.- Figure 10 
shows the variation of Strouhal number for different Richardson numbers. At 
the first place, the shedding frequency is observed to decrease with increased 
cooling of the cylinder below the fluid temperature. The fluid in the wake 
becomes denser due to cooling and therefore the vortices stay attached to the 
cylind(!r for long. 'I'lie hr()a<lening of tin; wake also helps in accomodating large 
size vortices which are shed in the downstrem at greater time intervals. For 
Richardson number upto 0.1 the shedding frequency steadily increses but for 
Ri > 0.15 the shedding frequency suddenly becomes zero i.e, the shedding is 
stopped. 


St 
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The temporal mean of the Nusselt number averaged over the entire cylinder 
surface has been plotted against the Richardson number in Fig.ll. The be- 
haviour of the average Nusselt number characteristic can also be explained on 
the basis of the above arguments. Over the range of positive Richardson num- 
ber, increasing magnitude of Ri tends to increase the average Nusselt number 
but for negative Ri no significant change in Nusselt number is observed because 
of longer residence time of the large size vortices at the back of the cylinder. 
The reduced shedding frequency in the regime of negative Richardson number 
tends to reduce the heat transfer. However beyond the critical Ri of 0.15 the 
vortices becomes small in size and thus the zone of poor heat transfer is re- 
duced. The increase in Nusselt number is totally attributed to the increase in 
the magnitude of Ri. 

The time average value of the local Nusselt number Nu^ has been plotted 
against the angle (j) in Fig.l2. The situation <^ = 0 corresponds to the stag- 
nation point facing the incoming fluid. The angle increases downstream till 
the tail of the cylinder is reached. The initial increase in local Nusselt number 
(upto (^=50) is due to the increase in velocity of the fluid. As the fluid moves 
downstream from (/> = 0 the flow area reduces and therefore the velocity rises. 
The Nu^ decreases with further increase in ^ due to the growth of boundary 
layer on the cylinder. The minimum point is the point of separation. It can 
be readily observed that the separation point is delayed with heating of the 
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cylinder. l‘'or lu'gat.ive Hi the h<;at. transfer increases rapidly after the point 
of separation bnt for positive Ri there is no significant rise in local heat trans- 
fer. In case of heating of the cylinder the buoyancy aids the forced convection 
and therefore the back-flow is greatly reduced. This reduces the mass of fluid 
present near the cylinder after the point of separation. There is therefore no 
significant increa.se in heat transfer. In case of cooling of the cylinder, due to 
the buoyancy the fluid returns back after separation and is responsible for the 
rapid increase in heat i.ransfer aft<ir the point of separation for negative Ri case. 

For comparison with earlier studies, the channel width has been enhanced. 
The blockage ratiod (D/II) of 0.1 and the far field boundary condition for the 
confining walls have been deployed. The Reynolds number of interest is 40 and 
the Richardson numbers are zero and unity. The comparisons of local Nusselt 
number distribution has been shown in Fig. 13(a). The results are in good 
agreement with the local Nusselt number distribution obtained by Badr [7] 
and Dennis etal [19]. Variation of Strouhal numbers as a function of Reynolds 
number has been shown in Fig. 13(b). All the computations have simulated 
flow past a cylinder in an infinite medium. The Strouhal numbers obtained in 
the present computation agree well with those due to Williamson [20] and Sa 
and Chang [21]. 







Chapter 5 

CONCLUSION 


A novel finite; voIuiik; algorithm has b<!on successfully extended to simulate a 
complex mixexl convection flow in a vertical channel with a built-in circular 
cylinder. 

T wo different flow regimes, depending on whether Richardson number is less 
than or greater than a critical value of 0.15 have been identified. The regime 
with Hi < 0.15 is characterized by broadening of the wake, advanced flow sep- 
aration and strong vortex shedding. The other regime {Ri > 0.15) represents 
increased separation delay and a zone of twin vortices attached to the cylinder. 
The Strouhal number reduces as the cylinder is cooled below the fluid tem- 
perature. Formation of large size vortices and reduced shedding frequency do 
not enhance heat transfer significantly for larger value of negative Richardson 
number. For Ri > 0.15 no vortex shedding takes place and the wake is also 
reduced in size. In this regime there is significant increase in average Nusselt 
number with increase in the magnitude of Richardson number. 
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